• Еще раз про множественные сравнения
  • Пример: Поведение песчанок в тесте открытое поле
  • perMANOVA
  • Условия применимости perMANOVA
  • Более подробная интерпретация результатов perMANOVA
  • Более сложные дизайны в perMANOVA

Вы сможете

  • Вспомнить основные дизайны для тестирования гипотез в рамках дисперсионного анализа
  • Применить функции R для тестирования гипотез с помощью perMANOVA

Еще раз про множественные сравнения

В чем опасность множественных сравнений?

\(\alpha\) — это вероятность совершить ошибку первого рода при тестировании гипотезы (= вероятность отвергнуть истинную нулевую гипотезу, = вероятность найти различия там, где их нет).

Обычно принимается, что \(H_0\) отвергают на уровне значимости \(\alpha = 0.05\).

Когда у нас два средних — все просто, сравнение всего одно. Естественно, вероятность совершить ошибку I рода для группы сравнений \(\alpha _{familywise}\) равна уровню значимости для единственного сравнения \(\alpha _{per\,comparison}\).

\(\alpha _{familywise} = \alpha _{per\,comparison}\)

Но если сравнений много, то растет вероятность совершить хотя бы одну ошибку I рода (найти различия там, где их нет).

Если сравнений много…

Например, если мы хотим попарно сравнить три значения, нам понадобится сделать 3 сравнения.

Пусть мы решили, что в каждом из сравнений будем использовать уровень значимости \(0.05\).

Тогда в каждом из сравнений вероятность совершить ошибку первого рода будет \(\alpha_{per\,comparison} = 0.05\).

Если сделать \(N\) тестов, то вероятность совершить хотя бы одну ошибку I рода в группе тестов (family-wise error rate, FWER) значительно возрастает.

\[\alpha_{familywise} = 1 - (1 - \alpha_{per\,comparison})^N\]

Чем больше сравнений, тем больше вероятность обнаружить различия там, где их на самом деле нет.

\[\alpha_{familywise} = 1 - (1 - \alpha_{per\,comparison})^N\]

В таблице даны значения \(\alpha _{familywise}\) для разного числа сравнений, если \(\alpha _{per\,comparison} = 0.05\):

Число средних Число сравнений αfamilywise
2 1 0.050
3 3 0.143
4 6 0.265
5 10 0.401

Для решения проблемы есть два подхода

  1. Взять более жесткий порог уровня значимости. (Например, ввести поправку Бонферрони, поправку Хольма-Бонферрони, воспользоваться процедурой Бенъямини-Хохберга).
  2. Изменить схему тестирования гипотезы — тестировать не три независимых гипотезы, а одну сложную (так это, например, происходит в ANOVA).

Поправка Бонферрони

\(\alpha^*_{per\,comparison} = \alpha_{familywise} /{N}\)

Это жесткий способ, т.к с возрастанием числа сравнений резко снижается уровень значимости и мощность теста. В результате растет риск не найти различий, где они есть.

Ниже \(\alpha _{per\,comparison}\) после поправки Бонферрони, сохраняющие \(\alpha _{familywise} = 0.05\):

Число средних Число сравнений αper comparison
2 1 0.050
3 3 0.017
4 6 0.008
5 10 0.005
6 15 0.003

Метод Хольма-Бонферрони

Метод Хольма-Бонферрони — это пошаговая процедура.

Чтобы зафиксировать \(FWER \le \alpha_{familywise}\):

  1. Сортируем в порядке возрастания \(N\) значений \(p\), полученные в тестах, и присваиваем каждой ранг \(j\) от 1 до \(N\):
    \[p_{1} \le p_{2} \le \cdots \le p_{N - 1} \le p_{N}\]
  2. Вводим поправку для значений \(p\)
    \[{p^*_{j}} = min{\{(N - j + 1) \cdot p_{j},\;1\}}\]
  3. Сравниваем с \(\alpha\). Отвергаем \(H_0\), если \(p^*_{j} < \alpha_{familywise}\)

Пример поправки Хольма-Бонферрони

В таблице ниже даны результаты нескольких сравнений (\(N = 5\)). С помощью поправки Хольма-Бонферрони для каждого \(p_j\) мы получим свое скорректированное значение \(p^*_{j}\).

Ранг(\(j\)) \(\mathbf{p_{j}}\) \(N - j + 1\) \(\mathbf{p^*_{j}} = min{\{(N - j + 1) \cdot p_{j},\;1\}}\) Отвергаем \(H_0\)?
1 0.001 5 0.005 Да
2 0.010 4 0.040 Да
3 0.035 3 0.105 Нет, и дальше везде сохраняем \(H_0\)
4 0.040 2 0.080
5 0.046 1 0.046

Дисперсионный анализ (повторение)

Классический дисперсионный анализ

Ronald Fisher

Рональд Элмер Фишер

Пусть имеется несколько градаций фактора A
(например, \(A _{1...3}\))

  • Почему появляется межгрупповая изменчивость, то есть разные средние значения для групп по фактору А?
  • Почему появляется внутригрупповая изменчивость, то есть разные значения y в группах?
A1 A2 A3
\(y_{11}, y_{12}, y_{13}\) \(y_{21}, y_{22}, y_{23}\) \(y_{31}, y_{32}, y_{33}\)
\(\bar{y}_{A_1}\) \(\bar{y}_{A_2}\) \(\bar{y}_{A_3}\) \(\bar{Y}_{gen}\)

Суммарная дисперсия может быть разложена на две составляющие

\(SS_{total} = SS_{between} + SS_{within}\)

Для тестирования гипотезы о влиянии фактора надо сравнить межгрупповую изменчивость (\(SS_{between}\)) и внутригрупповую (\(SS_{within}\)).

Для этого сравнения Фишер предложил статистику \(F = \frac{SS_{between} / (a - 1)}{SS_{within} / (N - a)}\), где \(a\) — число групп.

F-распределение

Если межгрупповая изменчивость равна внутригрупповой, то \(F\) принадлежит \(F\)-распределению с двумя параметрами \(df_{between} = a - 1\) и \(df_{within} = N - a\), где \(a\) — число классов, \(N\) — общее количество объектов в анализе.

Однофакторный дизайн

Выявляется влияние фактора А.

Многофакторный ортогональный дизайн

Выявляется влияние фактора А, B и их взаимодействия A\(\times\)B.

Иерархический дизайн

Некоторые дискретные предикторы могут быть иерархически соподчинены.

Выявляется влияние вложенного фактора B внутри градаций фактора A.

Влияние вложенного фактора B можно оценить, но чаще всего оно не представляет интереса для исследователя.

Рандомизированный полный блочный дизайн

Влияние блокриующего фактора тоже можно оценить, но часто оно не представляет интереса для исследователя. Обычно блокирующий фактор рассматривается как случайный.

Пример: Поведение песчанок в тесте открытое поле

Пример: Поведение песчанок в тесте открытое поле

Гипотеза: Разные виды песчанок различаются по поведению в тесте “Открытое поле”

Карликовая песчанка Монгольская песчанка Жирнохвостая песчанка

Виды:

  • Карликовая песчанка (Gerbillus gerbillus)
  • Монгольская песчанка (Meriones unguiculatus)
  • Жирнохвостая песчанка (Pachyuromys duprasi)





Фото: (1) XV8-Crisis; (2) Alastair Rae; (3) P.H.J. (Peter) Maas / www.petermaas.nl

Тест “открытое поле”

Пример: Поведение песчанок в тесте открытое поле

Гипотеза: Разные виды песчанок различаются по поведению в тесте “Открытое поле”

Карликовая песчанка Монгольская песчанка Жирнохвостая песчанка

Оценка поведения песчанок трех видов по семи признакам:

  • Время до выхода в квадрат открытого поля
  • Количество актов мочеиспускания
  • Количество актов дефекации
  • Количество пересеченных квадратов
  • Число вертикальных стоек
  • Количество актов смещенной активности
  • Время проведенное в центре квадрата открытого поля

Фото: (1) XV8-Crisis; (2) Alastair Rae; (3) P.H.J. (Peter) Maas / www.petermaas.nl

Данные наблюдений

pesch <- read.csv("data/pesch.csv", header = TRUE, sep = ";")
head(pesch)

Задание

  • Какую меру различия можно использовать с этими данными?
  • Как можно преобразовать данные?
  • Постройте ординацию объектов в осях MDS и раскрасьте точки в соответствии с видами

Решение

Поскольку измеренные признаки варьируют в разных масштабах, целесообразно логарифмировать данные.

options(digits = 4)
log_pesch <- pesch
log_pesch[, 3:ncol(pesch)] <- log(pesch[, 3:ncol(pesch)] + 1)
head(log_pesch)
##   Gender Species Time_to_entrance Urination Defecation Quadr_Number
## 1      f    karl            0.000    0.6931      0.000        3.871
## 2      f    karl            3.045    0.0000      1.386        5.762
## 3      m    karl            5.204    0.0000      0.000        5.182
## 4      f    karl            0.000    0.0000      0.000        3.497
## 5      f    karl            4.942    0.0000      0.000        5.328
## 6      m    karl            0.000    0.0000      1.099        3.664
##   Vert_Number Displ_Act Time_in_Centre
## 1       2.485     1.609          0.000
## 2       4.078     1.386          1.946
## 3       4.025     1.946          1.099
## 4       1.792     1.386          0.000
## 5       3.401     2.398          1.386
## 6       2.303     2.197          0.000

Решение

library(ggplot2)
library(vegan)
theme_set(theme_bw())

mds_pesch <- metaMDS(log_pesch[, 3:ncol(pesch)], distance = "euclidean")
mds_pesch <- as.data.frame(mds_pesch$points)
mds_pesch$Species <- pesch$Species

ggplot(mds_pesch, aes(x = MDS1, y = MDS2, colour = Species)) + 
  geom_point(size = 5)

Различаются ли виды песчанок по поведению?

Каким способом можно ответить на этот вопрос?

  • Мы могли бы взять каждый из поведенческих признаков отдельно и провести одномерный однофакторный дисперсионный анализ.
  • Но нас интересует поведение в целом — нужен многомерный анализ.

Методы выявления различий между группами для многомерных данных

ANOVA разработан для одномерных данных

Что делать если мы хотим оценивать объект по многим признакам сразу?

Примеры:

  • Сообщество как целое (много видов)
  • Поведение как целое (много элементов)
  • Метаболом как целое (много метаболитов)
  • Ответы респондентов на множество вопросов в анкетах

Варианты решений:

  • MANOVA (Fisher, 1925, Wilks, 1932)
  • distance-based Redundancy Analysis (db-RDA) (Legendre, Anderson, 1999)
  • perMANOVA (Anderson, 2001; McArdle, Anderson, 2001)

Многомерный дисперсионный анализ (MANOVA)

Давно разработан параметрический метод MANOVA (Multivariate Analysis Of Variance). Он дает возможность проводить анализ аналогичный ANOVA. В основе MANOVA лежат представление о многомерном нормальном распределении и расстояниях между центроидами.

В MANOVA сравниваются:

  • отклонения точек от групповых центроидов (аналог \(SS_{within}\))
  • отклонения групповых центроидов от общего центроида (аналог \(SS_{between}\)).

MANOVA

Anderson, 2001





Ограничения MANOVA:

  • Многомерная нормальность распределения
  • Гомогенность дисперсий

Permutational Multivariate Analysis of Variance (perMANOVA)

Marti Anderson

Марти Джейн Андерсон





Anderson 2001

Теорема

Сумма квадратов отклонений от центроидов равна сумме квадратов взаимных расстояний, деленной на число объектов

Для Евклидовых расстояний эта закономерность была известна давно (например, Kendall, Stuart 1963).

Centroid and distanceAnderson, 2001





В случае Евклидова расстояния (именно его имплицитно использует MANOVA) центроиды найти очень просто — это средние значения соответствующих координат. Поэтому обычно сначала непосредственно вычисляли центроиды, и затем — сумму квадратов отклонений от них.

Для других мер различия центроиды найти гораздо сложнее. Например, для коэффициента Брея-Куртиса (не метрика), среднее значение не будет соответствовать центроиду.

Марти Андерсон показала, что можно обойтись без вычисления центроидов и для других мер различия

Centroid and distance

Anderson, 2001

MANOVA (Евклидово расстояние)

MANOVA

Anderson, 2001

perMANOVA (любой коэффициент)

perMANOVA


Можно непосредственно из матрицы любых коэффициентов различия найти и общую и внутригрупповые суммы квадратов

 Anderson, 2001

Разложение дисперсии становится очень простым

distances-total-within-between Anderson, 2001



Пусть всего N элементов, принадлежащих a группам по n элементов в каждой, d - расстояние между i-тым и j-тым объектами, \(\epsilon\) - 1 если объекты i и j из одной группы и 0, если из разных.

\(SS_{total} = \frac{1}{N}\sum \sum {d_{ij}^2}\)
Сумма квадратов взаимных расстояний — это сумма квадратов субдиагональных элементов, деленная на число объектов N.

\(SS_{within} = \frac{1}{n}\sum \sum {d_{ij}^2 \cdot \epsilon_{ij}}\)
Внутригрупповая сумма квадратов — это сумма всех сумм квадратов расстояний между элементами для каждой группы, деленная на n число объектов в группе.

Тогда межгрупповая сумма квадратов \(SS_{between} = SS_{total} - SS_{within}\)

Псевдо-F статистика

\[F = \frac{SS_{between} / (a - 1)}{SS_{within}/(N - a)}\]

Для оценки значимости псевдо-F используется пермутационная процедура:

  • Случайным образом перетасовываются строки исходной матрицы
  • После каждого акта пермутации вычисляется \(F_{perm}\)
  • Уровень значимости
  • Внимание! Для иерархического дизайна процедура пермутации имеет свои особенности (обсудим позднее).

\[p = \frac {N_{F_{perm} \ge F}} {N_{permutations}}\]

perMANOVA в примере

Применим метод perMANOVA (функция adonis())

library(vegan)
permanova_pesch <- adonis(log_pesch[3:9] ~ Species, data = log_pesch,
                          method = "euclidean")
permanova_pesch
## 
## Call:
## adonis(formula = log_pesch[3:9] ~ Species, data = log_pesch,      method = "euclidean") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)  
## Species    2        60   30.02    3.58 0.237  0.012 *
## Residuals 23       193    8.39         0.763         
## Total     25       253                 1.000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Мы видим традиционную для ANOVA таблицу результатов. Что здесь что?

Результаты perMANOVA

## 
## Call:
## adonis(formula = log_pesch[3:9] ~ Species, data = log_pesch,      method = "euclidean") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)  
## Species    2        60   30.02    3.58 0.237  0.012 *
## Residuals 23       193    8.39         0.763         
## Total     25       253                 1.000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Поведение разных видов песчанок в тесте открытое поле значимо различалось (perMANOVA, p = 0.012).

Условия применимости perMANOVA

Условия применимости perMANOVA

  1. Требуется равенство внутригрупповых дисперсий (гомоскедастичность).

  2. Желательно использование сбалансированных данных (с равным объемом групп), т.к. этом случае perMANOVA устойчив к гетерогенности дисперсии (Anderson, Walsh, 2013).

Проверка равенства внутригрупповых дисперсий

Для проверки можно использовать функцию betadisper(), изначально предназначенную для сравнения \(\beta\)-разнообразия сообществ.

Функция betadisper() вычисляет внутригрупповые центроиды и координаты точек в пространстве главных координат (Principal Coordinate Analysis = PCoA = metric MDS).

Затем, значимость различий отклонений от центроидов в разных группах проверяется с помощью процедуры PERMDISP2.

dist_pesch <- vegdist(log_pesch[,3:ncol(pesch)], method  = "euclidean")
PCO_pesch <- betadisper(dist_pesch, log_pesch$Species)

График ординации PCoA

Объект, возвращаемый betadisper(), позволяет также нарисовать наши объекты в пространстве главных координат (PCoA).

plot(PCO_pesch, main = "PCoA ordination")

Процедура PERMDISP2 для проверки равенства внутригрупповых дисперсий

Процедура PERMDISP2 реализована в пакете vegan в функции anova.betadisper(). Это многомерный аналог теста Левина на гомогенность дисперсий в группах, который иногда используется для проверки условий применимости дисперсионного анализа.

anova(PCO_pesch)
## Analysis of Variance Table
## 
## Response: Distances
##           Df Sum Sq Mean Sq F value Pr(>F)
## Groups     2    9.7    4.86    2.05   0.15
## Residuals 23   54.5    2.37
  • Не выявлено значимых различий разброса внутригрупповых расстояний.

Для визуализации можно нарисовать боксплот

boxplot(PCO_pesch)

Более подробная интерпретация результатов perMANOVA

Post hoc тесты в perMANOVA

На приведенной ординации видно, что точки, соответствующие Монгольским песчанкам расположены отдельно от остальных.

Для выявления попарных различий нужны попарные сравнения.

Попарные сравнения как замена post hoc тесту

Внимание! В пакете vegan пост хок тест не реализован. Но мы можем сделать простейшую его версию самостоятельно.

Карликовая песчанка Монгольская песчанка Жирнохвостая песчанка


Проведем попарные сравнения между группами, то есть:

  • Карликовые vs. Монгольские
  • Карликовые vs. Жирнохвостые
  • Монгольские vs. Жирнохвостые




Фото: (1) XV8-Crisis; (2) Alastair Rae; (3) P.H.J. (Peter) Maas / www.petermaas.nl

Функция для попарных сравнений perMANOVA

pairwise_permanova <- function(dat, group, strata = NULL, ...){
  pair <- combn(unique(as.character(group)), 2)
  ncomb <- ncol(pair)
  res <- rep(NA, ncomb)
  for (i in 1:ncomb) {
    filter <- group %in% pair[, i]
    if(is.null(strata)){
      posthoc <- adonis(dat[filter, ] ~ group[filter], ...)$aov.tab$Pr[1]
    } else {
      posthoc <- adonis(dat[filter, ] ~ group[filter], 
                        strata = strata[filter], ...)$aov.tab$Pr[1]
    }
    res[i] <- posthoc
    names(res)[i] <- paste(pair[, i], collapse = " vs. ")
  }
  return(res)
}

Результаты попарных сравнений

Карликовая песчанка Монгольская песчанка Жирнохвостая песчанка
p_vals <- pairwise_permanova(
  dat = log_pesch[, -c(1:2)], group = log_pesch$Species, 
  method = "euclidean", permutations=9999)
p_vals
##         karl vs. mongol   karl vs. zhirnokhvost mongol vs. zhirnokhvost 
##                  0.0011                  0.4029                  0.0031

Это все? Пишем статью?



Фото: (1) XV8-Crisis; (2) Alastair Rae; (3) P.H.J. (Peter) Maas / www.petermaas.nl

Поправка на множественные сравнения

Мы делали три пары сравнений — нужно ввести поправку для уровня значимости. Будем считать значимыми различия в тех сравнениях, где после введения поправки \(p < 0.05\).

Можно посчитать скорректированные значения уровня значимости \(p\) с учетом поправки Хольма-Бонферрони.

p.adjust(p_vals, method = "holm")
##         karl vs. mongol   karl vs. zhirnokhvost mongol vs. zhirnokhvost 
##                  0.0033                  0.4029                  0.0062

Более сложные дизайны в perMANOVA

Многофакторный дизайн в perMANOVA

Выясним, влияет ли пол и вид песчанок на поведение.

Отфильтруем исходные данные (в случае с жирнохвостыми песчанками были изучены только самки)

log_pesch2 <- log_pesch[log_pesch$Species != "zhirnokhvost", ]

Проведем двухфакторный анализ perMANOVA

Различается ли поведение песчанок в зависимости от видовой принадлежности и пола?

twofact_pesch <- adonis(log_pesch2[,3:ncol(pesch)] ~ Gender * Species,
                        data = log_pesch2, method = "euclidean")
twofact_pesch
## 
## Call:
## adonis(formula = log_pesch2[, 3:ncol(pesch)] ~ Gender * Species,      data = log_pesch2, method = "euclidean") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)   
## Gender          1       7.2     7.2    1.13 0.049  0.292   
## Species         1      45.5    45.5    7.16 0.311  0.004 **
## Gender:Species  1       4.5     4.5    0.72 0.031  0.529   
## Residuals      14      88.9     6.4         0.608          
## Total          17     146.1                 1.000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Здесь возможен иерархический дизайн

Различается ли поведение самцов и самок у этих видов песчанок?

nested_pesch <- adonis(log_pesch2[, 3:ncol(pesch)] ~ Gender, 
                       data = log_pesch2, strata = log_pesch2$Species, 
                       method = "euclidean")
nested_pesch
## 
## Call:
## adonis(formula = log_pesch2[, 3:ncol(pesch)] ~ Gender, data = log_pesch2,      method = "euclidean", strata = log_pesch2$Species) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)
## Gender     1       7.2    7.21    0.83 0.049    0.3
## Residuals 16     138.9    8.68         0.951       
## Total     17     146.1                 1.000

Внимание! Пермутации производятся только в пределах группирующего фактора (аргумент strata)

Задание

  • Создайте датафрейм из файла simulated_data.csv (Это данные симулированные по алгоритму, приведенному в справке по функции adonis())
  • В этом датафрейме записано обилие двух видов на экспериментальных площадках двух типов: без добавления и с добавлением NO3, по 6 повторностей в каждом эксперименте.
    Эксперименты были независимо проведены на 3 полях.
  • Оцените, зависит ли структура совместного поселения этих двух видов от концентрации NO3.

Решение

com <- read.csv("data/simulated_data.csv", sep = ',', header = T)

# Ошибочный дизай
com_permanova <- adonis(com[,1:2] ~ com$NO3)

# Правильный дизайн
com_permanova2 <- adonis(com[,1:2] ~ com$NO3, strata = com$field)

Summary

  • При одновременном тестировании нескольких гипотез растет вероятность ошибки I рода.
  • Чтобы контролировать вероятность ошибки I рода, либо используют более жесткий уровень значимости, либо тестируют сложную гипотезу вместо нескольких простых.
  • perMANOVA дает возможность тестировать сложные гипотезы в отношении явлений, описанных по многим переменным (т.е. на многомерных данных).
  • В perMANOVA можно использовать любые коэффициенты различия.
  • Для применения perMANOVA требуется равенство разбросов точек между центроидами их групп, но при равных объемах групп анализ устойчив к отклонениям от этого условия.
  • При использовании perMANOVA важно не запутаться в дизайне.

Другие программы

  • Primer 6.0 + PERMANOVA Коммерческий продукт.
  • PAST Некоммерческая программа. Здесь метод называется NPMANOVA.
  • Оригинальная программа М. Андерсон PERMANOVA и PERMDISP.

Что почитать

  • Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology, 26: 32–46.
  • Anderson, M.J. 2005. PERMANOVA: a FORTRAN computer program for permutational multivariate analysis of variance. Department of Statistics, University of Auckland, New Zealand.
  • Anderson, M.J. (2004). PERMDISP: a FORTRAN computer program for permutational analysis of multivariate dispersions (for any two-factor ANOVA design) using permutation tests. Department of Statistics, University of Auckland, New Zealand.
  • Legendre P., Legendre L. (2012) Numerical ecology. Second english edition. Elsevier, Amsterdam.